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ABSTRACT 


This  project,  supported  by  a  grant  from  the  Air  Force  Office  of  Scientific  Research,  was 
primarily  aimed  at  the  design  of  novel  algorithms  for  the  restoration  and  super-resolution 
processing  of  imagery  data  to  improve  the  resolution  in  images  acquired  from  practical  sensing 
operations.  Due  to  the  diffraction  limits  underlying  these  operations,  the  image  recorded  at  the 
output  of  the  imaging  system  is  a  low-pass  filtered  version  of  the  original  scene  and  hence 
recovery  of  the  lost  information  contributing  to  the  finer  details  is  required  to  produce  significant 
resolution  enhancement.  While  traditional  image  restoration  procedures  (based  on  deconvolution 
and  inverse  filtering  approaches)  attempt  mainly  reconstruction  of  the  passband,  super-resolution 
algorithms  attempt  to  provide  not  only  passband  restoration  but  also  some  degree  of  spectral 
extrapolation  thus  enabling  to  restore  the  high  frequency  spatial  amplitude  variations  relating  to 
the  spatial  resolution  of  the  sensor  and  lost  due  to  diffraction-limited  imaging.  These  algorithms 
are  typically  iterative  in  nature  and  implement  nonlinear  signal  processing  operations.  Two 
distinct  approaches  that  have  resulted  in  powerful  super-resolution  algorithms  are  based  on 
statistical  optimization  arguments  and  set-theoretic  estimation  procedures.  The  principal 
objectives  in  this  project  were  to  develop  and  evaluate  specific  techniques  for  developing  new 
processing  algorithms  by  combining  the  strong  points  of  the  two  approaches.  The  principal 
outcomes  from  this  work  include  the  following: 

•  Systematic  procedures  for  extracting  and  modeling  scene-derived  information  sets  for 
projection-based  set-theoretic  super-resolution  processing; 

•  Design  of  hybrid  processing  algorithms  that  integrate  Projection  Onto  Convex  Sets 
(POCS)  iterations  with  Maximum  Likelihood  (ML)  estimation  procedures  to  yield 
superior  restoration  and  super-resolution  performance. 

In  order  to  increase  the  relevance  of  the  work  conducted  in  this  project  to  the  U.S'.  Air  Force 
needs,  a  close  liaison  was  maintained  with  the  Air  Force  Wright  Laboratory  Armament 
Directorate  (at  Eglin  AFB).  For  proof  of  concept  demonstrations,  the  tailoring  of  the  algorithms 
as  well  as  the  processing  studies  were  conducted  with  a  special  focus  on  achieving  resolution 
enhancements  in  the  passive  millimeter- wave  (PMMW)  images  acquired  by  the  Air  Force 
Wright  Laboratory  Armament  Directorate  with  the  state-of-the-art  radiometers  being  built  by 
them. 
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1.  INTRODUCTION 


1.1  Principal  Objectives  of  the  Project 

Imagery  data  acquired  from  practical  sensing  operations  typically  suffer  from  poor  spatial 
resolution  due  to  the  finite  size  of  antenna  or  lens  that  makes  up  the  imaging  system  and  the 
consequent  imposition  of  the  underlying  diffraction  limits.  Hence,  some  form  of  post-detection 
image  processing  directed  to  enhancing  the  resolution  will  often  be  needed  before  the  obtained 
measurements  can  be  used  reliably  for  any  image  exploitation  and  decision-making  objectives. 
Image  restoration  has  a  considerably  long  history  and  several  powerful  procedures  have  been 
developed  to  restore  degraded  images  [1].  However,  since  the  image  recorded  at  the  output  of  the 
imaging  system  is  a  low-pass  filtered  version  of  the  original  scene,  recovery  of  the  lost 
information  contributing  to  the  finer  details  is  indeed  necessary  to  produce  significant  resolution 
enhancement.  Recognizing  this  fact,  researchers  have  directed  their  attention  to  the  design  of 
efficient  processing  schemes,  commonly  referred  to  as  super-resolution  algorithms,  that  achieve 
not  only  a  satisfactory  restoration  of  the  passband  (by  attempting  to  reverse  the  effects  of 
convolution  with  the  sensor  point  spread  function  within  the  image  bandwidth)  but  also  some 
degree  of  spectral  extrapolation  to  recover  the  lost  frequencies  beyond  the  diffraction-limited 
cutoff. 

The  primary  goal  in  this  project  was  the  development  of  novel  approaches  and 
procedures  for  the  restoration  and  super-resolution  of  diffraction-limited  imagery  data  in  order  to 
achieve  resolution  enhancements  that  facilitate  improved  surveillance  and  tracking  of  targets  of 
interest.  In  practical  applications,  potential  advantages  that  could  be  achieved  from  such 
resolution  enhancements  beyond  those  dictated  from  practical  limits  (such  as  the  one  given  by 
the  Rayleigh  criterion  [2])  are: 
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(i)  benefits  in  range  of  measurement  -  scene  details  that  are  potentially  available  from  a  closer 
range  can  be  obtained  from  longer  range  measurements; 

(ii)  benefits  in  sensor  dimensions  -  resolution  levels  that  are  potentially  obtainable  from  a  larger 
aperture  or  a  bigger  antenna  can  be  achieved  from  less  expensive  or  smaller  size  sensors  (  a 
typical  payoff  being  to  gain  the  effect  of  deploying  a  big  antenna  on  a  small  airframe). 

Our  investigations  towards  this  goal  hence  focused  on  several  specific  topics  that  included 
the  following: 

•  development  of  novel  approaches  based  on  statistical  modeling  and  convex  constraint  sets  for 
resolution  enhancement  in  the  presence  of  noise  and  clutter; 

•  mathematical  modeling  of  constraint  sets  and  projection  operators  for  specific  a  priori 
information  about  the  object  or  scene  being  imaged  in  order  to  facilitate  inclusion  in  the 
design  of  systematic  restoration  and  super-resolution  procedures  based  on  the  Projection- 
Onto-Convex-Sets  (POCS)  approach; 

•  investigation  of  methods  for  combining  POCS-based  approaches  with  maximum  likelihood 
(ML)  restoration  methods  for  designing  hybrid  POCS-Bayesian  super-resolution  procedures; 

Due  to  the  fact  that  investigations  on  resolution  enhancement  have  both  military  and  non¬ 
military  applications,  the  studies  were  conducted  using  a  general  framework  and  employing 
mathematical  models.  However,  in  order  to  provide  validation  of  the  studies  conducted  and  the 
results  obtained  to  military  environments,  the  work  was  performed  in  close  collaboration  with 
the  Air  Force  Wright  Laboratory  Armament  Directorate  personnel.  In  particular,  the  various 
simulation  experiments  conducted  to  demonstrate  the  efficiency  of  the  algorithms  developed  in 
this  project  were  appropriately  tailored  to  process  passive  millimeter-wave  (PMMW)  imagery 
data  acquired  from  state-of-the-art  PMMW  radiometers  being  built  by  the  Advanced  Sensors 
Group  at  the  Air  Force  Wright  Laboratory  Armament  Directorate. 
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1.2  Restoration  and  Super-resolution  Processing  of  Sensor  Acquired  Data 

Imaging  sensors,  no  matter  how  carefully  designed,  have  certain  limitations  arising 
mainly  from  physical  constraints,  the  principal  one  being  poor  resolution  in  the  acquired 
imagery.  The  diffraction-limited  angular  resolution,  Q ,  of  an  incoherent  imaging  system  [2]  is 
given  by 

A 

e-Mi- 

where  A  is  the  effective  wavelength  of  imaging  and  D  is  the  diameter  of  the  limiting  aperture  of 
the  antenna  or  lens.  As  A  increases,  the  achievable  angular  resolution  decreases,  i.e.  6 ,  the  size 
of  the  angle  between  two  resolvable  points,  increases.  It  may  be  noted  that  the  wavelength  of  a 
synthetic  aperture  radar  (SAR)  operating  at  1  GHz  is  about  1  inch  long  and  one  needs  an  antenna 
as  big  as  40  ft  wide  in  order  to  achieve  a  resolution  requirement  of  being  able  to  able  to 
distinguish  points  in  a  scene  separated  by  about  1  meter  at  a  distance  of  1  Km.  For  a  typical 
passive  millimeterwave  (PMMW)  sensor  with  a  1  ft  diameter  antenna  and  operating  at  95  GHz, 
the  angular  resolution  is  only  about  1 0  mrad,  which  translates  into  a  spatial  resolution  of  about 
10  meters  at  a  distance  of  1  Km.  Some  recent  studies  have  also  established  that  for  ensuring 
reasonably  adequate  angular  resolution  (typically  of  the  order  of  4  mrad),  a  95  GHz  PMMW 
imaging  system  with  a  sensor  depression  angle  of  60°  -  80°  needs  to  be  confined  to  very  low 
operational  altitudes  (of  the  order  of  75-100  meters),  which  puts  inordinate  demands  on  the 
surveillance  and  guidance  schemes  to  facilitate  such  requirements.  Similar  resolution  limitations 
and  the  consequent  requirements  on  operational  conditions  (some  of  which  may  be  clearly 
impossible  to  satisfy  for  tactical  missions  with  reliability  and  survivability  constraints)  exist  for 
other  types  of  sensing  modalities  as  well.  This  in  turn  demands  use  of  clever  image  processing 
techniques  to  process  the  acquired  imagery  data  and  obtain  improved  resolution  in  the  processed 
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images. 


The  fundamental  problem  underlying  the  sensing  operation  is  the  "low-pass"  filtering 
effect  due  to  the  finite  size  of  the  antenna  or  lens  that  makes  up  the  imaging  system  and  the 
consequent  imposition  of  the  underlying  diffraction  limits.  Hence  the  image  recorded  at  the 
output  of  the  imaging  system  is  a  low-pass  filtered  version  of  the  original  scene.  The  portions  of 
the  scene  that  are  lost  by  the  imaging  system  are  the  fine  details  (high  frequency  spectral 
components)  that  accurately  describe  the  objects  in  the  scene,  which  also  are  critical  for  reliable 
detection  and  classification  of  targets  of  interest  in  the  scene.  Hence  some  form  of  image 
processing  to  restore  the  details  and  improve  the  resolution  of  the  image  will  invariably  be 
needed.  Traditional  image  restoration  procedures  (based  on  deconvolution  and  inverse  filtering 
approaches)  attempt  mainly  at  reconstruction  of  the  passband  and  possibly  elimination  of  effects 
of  additive  noise  components.  These  hence  have  only  limited  resolution  enhancement 
capabilities.  Greater  resolution  improvements  can  only  be  achieved  through  a  class  of  more 
sophisticated  algorithms,  called  super-resolution  algorithms,  which  provide  not  only  passband 
reconstruction  but  also  some  degree  of  spectra!  extrapolation,  thus  enabling  to  restore  the  high 
frequency  spatial  amplitude  variations  relating  to  the  spatial  resolution  of  the  sensor  and  lost 
through  the  filtering  effects  of  the  seeker  antenna  pattern.  A  tactful  utilization  of  the  imaging 
instrument's  characteristics  and  any  a  priori  knowledge  of  the  features  of  the  target  together  with 
an  appropriately  crafted  nonlinear  processing  scheme  is  what  gives  the  capability  to  these 
algorithms  for  super-resolving  the  input  image  by  extrapolating  beyond  the  passband  range  and 
thus  extending  the  image  bandwidth  beyond  the  diffraction  limit  of  the  imaging  sensor. 

For  application  in  surveillance  environments,  it  must  be  emphasized  that  super-resolution 
is  a  post-processing  operation  applied  to  the  acquired  imagery  and  consequently  is  much  less 


8 


expensive  compared  to  improving  the  imaging  system  for  desired  resolution.  As  an  example,  it 
may  be  noted  that  for  visual  imagery  acquired  from  space-borne  platforms,  some  studies  indicate 
that  the  cost  of  camera  payload  increases  as  the  inverse  2.5  power  of  the  resolution.  Hence  a 
possible  two-fold  improvement  in  resolution  by  super-resolution  processing  in  this  application 
roughly  translates  into  a  reduction  in  the  cost  of  the  sensor  by  more  than  5  times.  Similar 
relations  also  exist  for  sensors  operating  in  other  spectral  ranges  (due  to  the  relation  between 
resolution  and  antenna  size),  confirming  the  cost  effectiveness  of  employing  super-resolution 
algorithms.  The  principal  goal  of  super-resolution  processing  in  a  multispectral  surveillance 
environment  is  hence  to  obtain  an  image  of  a  target  of  interest  via  post-processing  that  is 
equivalent  to  one  acquired  through  a  more  expensive  larger  aperture  sensor. 

1.3  Basic  Considerations  in  Design  of  Super-resolution  Processing  Algorithms 

The  idea  of  recreating  the  spectral  components  that  are  removed  by  the  imaging  process 
and  hence  are  not  present  in  the  image  available  for  processing  may  pose  some  conceptual 
difficulties,  which  may  lead  one  to  suspect  whether  super-resolution  is  indeed  possible. 
Fortunately  there  exist  sound  mathematical  arguments  confirming  the  possibility  of  spectral 

■  i 

extrapolation.  The  primary  justification  comes  from  the  Analytic  Continuation  Theorem  and  the 
property  that  when  an  object  has  finite  spatial  extent  its  frequency  spectrum  is  analytic  [3].  Due 
to  the  property  that  a  finite  segment  of  any  analytic  function  in  principle  determines  the  whole 
function  uniquely,  it  can  be  readily  proved  that  knowledge  of  the  passband  spectrum  of  the 
object  allows  a  unique  continuation  of  the  spectrum  beyond  the  diffraction  limit  imposed  by  the 
imaging  system.  It  must  be  emphasized  that  the  limited  spatial  extent  of  the  object  is  critical  in 
providing  this  capability  for  extrapolation  in  the  frequency  domain. 
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Fundamental  to  the  reliable  estimation  of  the  high  frequencies  (that  underly  the  super¬ 
resolution  objectives)  is  the  utilization  of  a  priori  known  information  during  the  processing 
steps.  Indeed,  since  image  restoration  is  inherently  an  ill-conditioned  inverse  problem,  it  is  long 
realized  that  the  quality  of  restoration  and  the  extent  of  achievable  super-resolution  depend  on 
the  accuracy  and  the  amount  of  a  priori  information  the  processing  scheme  can  successlully 
employ.  Some  of  the  earliest  procedures  proposed  for  image  super-resolution  [4,5]  are  based  on 
intelligently  utilizing  the  known  information  in  the  form  of  constraints  during  the  iterative 
processing  steps.  In  practice,  one  finds  that  a  lot  of  information  typically  exists  to  assist  guiding 
the  restoration  process  in  a  satisfactory  direction.  This  includes,  (i)  Information  about  the 
solution ,  such  as  non-negativity  of  pixel  intensities,  spatial  extent  of  objects  present  in  the  scene, 
known  ranges  for  signal  intensities,  etc.,  (ii)  Information  about  the  imaging  system ,  such  as 
sensor  phenomenology,  aperture  size  and  the  resulting  shape  of  the  optical  transfer  function 
(OTF),  etc.,  and  (iii)  Information  about  the  imaging  conditions,  such  as  the  extent  of  observation 
noise,  object  motion  during  imaging,  etc.  In  addition,  as  will  be  shown  in  this  paper,  a  variety  of 
information  can  be  extracted  from  the  image  that  is  being  processed,  which  we  will  call  scepe- 
derived  information,  and  this  can  further  assist  in  the  restoration  operations. 

Does  availability  of  more  information  automatically  translate  into  better  restoration  and 
more  super-resolution?  Since  the  primary  objective  is  the  creation  of  data  not  present  in  the 
image  being  processed,  there  is  major  concern  on  the  fidelity  of  the  estimated  high  frequencies 
(one  may  note  that  addition  of  random  noise  can  also  lead  to  an  expansion  of  image  bandwidth). 
Utilization  of  more  and  diverse  types  of  information  consequently  provides  a  mechanism  for 
cross-validation  thus  resulting  in  a  faithful  re-creation  of  the  true  frequencies  present  in  the 
object  or  scene  being  restored.  Employing  all  of  the  available  information  does  not  come 
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without  a  price,  however.  More  information  used  during  the  iterative  processing  steps  almost 
always  translates  into  increased  computational  complexity,  which  in  turn  could  contribute  to 
poorer  convergence  speed  of  the  super-resolution  algorithm.  It  is  easy  to  see  that  if  two  different 
types  of  information  that  are  desired  to  be  utilized  are  not  appropriately  modeled  for  inclusion  in 
the  algorithm  (as  an  example,  two  information  sets  modeled  such  that  they  do  not  have  a  proper 
intersection),  it  will  have  detrimental  effects  on  the  algorithm  convergence.  Consequently, 
information  analysis  and  information  modeling  are  critical  steps  before  subjecting  the  image  to 
super-resolution  processing. 

Varying  by  the  extent  to  which  a  priori  knowledge  can  be  incorporated  in  algorithm 
development,  there  have  been  introduced  into  the  literature  a  large  number  of  image  restoration 
approaches  and  algorithms  too  vast  to  describe  or  reference  here.  One  may  refer  to  some  r.ecent 
survey  papers  [6,7]  for  a  review  of  the  extensive  activity  on  this  topic.  In  the  context  of  the 
objectives  of  the  present  project  however,  it  is  important  to  recognize  that  only  a  small  subset  of 
the  approaches  that  are  developed  for  image  restoration  have  received  some  interest  for  their 
super-resolution  capabilities,  i.e.  possible  spectrum  extrapolation  performance.  One  may  note 
that  not  all  image  restoration  methods  provide  the  capability  for  super-resolving.  In  fact,  a 
majority  of  existing  schemes  may  perform  decent  passband  restoration,  but  provide  no 
bandwidth  extension  at  all. 

The  various  approaches  in  general  attempt  to  code  the  a  priori  knowledge  to  be  used  by 
specifying  an  object  model  or  a  set  of  constraint  functions,  and  further  employ  an  appropriate 
optimization  criterion  to  guide  in  the  search  for  the  best  estimate  of  the  object.  A  convenient  way 
of  classifying  the  resulting  algorithms  is  into  iterative  and  noniterative  (or  direct)  schemes.  Fig. 
la  and  lb  depict  schematically  these  two  basic  approaches.  Noniterative  approaches  generally 
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attempt  to  implement  an  inverse  filtering  operation  (without  necessarily  performing  the 
computation  of  the  inverse  of  the  PSF  matrix  H,  however)  and  have  poor  noise  characteristics. 
All  required  computations  and  any  possible  use  of  constraint  functions  are  applied  in  one  step.  In 
contrast,  iterative  methods  apply  the  constraints  in  a  distributed  fashion  as  the  solution 
progresses  and  hence  the  computations  at  each  iteration  will  be  generally  less  intensive  than  the 
single-step  computation  of  noniterative  approaches.  Some  additional  advantages  of  iterative 
techniques  are  that,  (1)  they  are  more  robust  to  errors  in  the  modeling  of  the  image  formation 
process  (uncertainties  in  the  model  for  the  sensor  point  spread  function,  for  instance),  (2)  the 
solution  process  can  be  better  monitored  as  it  progresses,  (3)  constraints  can  be  utilized  to  better 
control  the  effects  of  noise  (and  possibly  clutter),  and  (4)  can  be  tailored  to  offset  sensor 
nonlinearities.  The  disadvantages  of  these  methods  generally  are,  (1)  increased  computation 
time,  and  (2)  need  for  proving  convergence  of  the  iterative  scheme.  Despite  these  disadvantages, 
iterative  methods  are  generally  the  preferred  approach  due  to  their  numerous  advantages  and  also 
since  the  iteration  can  be  terminated  once  a  solution  of  a  reasonable  quality  is  achieved. 
Consequently,  most  of  the  research  in  tailoring  super-resolution  algorithms  is  taking  place  at 
present  times  in  formulating  iterative  procedures. 


Fig.  la.  Schematic  of  Noniterative  ( Direct)  Super-resolution 
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A  priori  knowledge  and  constraints 
Fig.  lb.  Schematic  of Iterative  Super-resolution 
1.4  Approaches  and  Methods  for  Design  of  Super-resolution  Algorithms 

For  the  design  of  specific  digital  processing  algorithms  for  image  super-resolution,  a 
number  of  different  approaches  have  been  proposed  in  the  past,  which  have  lead  to  explicit 
procedures  for  both  iterative  and  non-iterative  constructions  of  the  restored  image  from  its 
degraded  version.  Perhaps  one  of  the  earliest  procedures,  suggested  by  Gerchberg  [4]  and 
Papoulis  [5],  involves  alternately  applying  constraints  in  the  space-domain  and  in  the  frequency- 
domain  in  the  quest  for  reconstructing  the  unknown  portion  of  the  frequency  spectrum  from  a 
limited  known  portion,  viz.  the  passband  of  the  sensor.  Much  of  the  activity  in  recent  times  has 
centered  on  statistical  optimization  methods.  In  this  approach  one  employs  an  iterative  scheme 
for  constructing  estimates  of  the  object  imaged,  either  in  the  form  of  a  maximum  likelihood 
(ML)  estimate  [8-10]  or  a  maximum  a  posteriori  (MAP)  estimate  [11-13],  from  using  the 
acquired  image  and  an  OTF  model  for  the  sensor.  A  specific  iterative  updating  scheme,  which  is 
a  special  form  of  implementing  the  well  known  Expectation-Maximization  (EM)  algorithm  and 
is  popularly  referred  to  as  Richardson- Lucy  iteration,  attempts  to  improve  the  value  of  the 
likelihood  function  during  successive  iterations.  This  scheme  has  received  considerable  attention 
in  the  past  and  several  studies  (too  numerous  to  mention  individually)  have  been  directed  to 
understanding  the  convergence  and  other  properties  of  the  iterative  process.  Several 
modifications  and  refinements  of  these  algorithms  have  also  been  developed  in  the  recent  past  to 
ensure  good  performance  even  under  non-ideal  conditions  such  as  imperfect  knowledge  of 
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sensor  OTF  and  high  levels  of  noise  corrupting  the  acquired  image  [14,15].  A  few  non-iterative 
algorithms  have  also  been  developed.  Notable  among  these  is  the  algorithm  due  to  Darling  et  al. 
[16],  which  performs  regularized  approximation  in  a  weighted  Hilbert  space  by  utilizing  a  priori 
information.  A  more  recent  effort  in  this  direction  is  due  to  Walsh  and  Nielsen-Delaney  [17]  who 
employed  a  method  based  on  singular  value  decomposition.  It  is  generally  agreed,  however,  that 
non-iterative  methods  do  not  offer  the  same  degree  of  flexibility  provided  by  iterative  methods  in 
utilizing  the  known  a  priori  information  during  processing,  which  is  at  the  heart  of  super¬ 
resolution  constructions.  Extensive  comparison  studies  that  document  the  performance  benefits 
resulting  from  the  use  of  iterative  algorithms  when  compared  with  other  approaches  for  the 
super-resolution  processing  of  specific  types  of  images  are  also  available  [18,35]. 

The  effectiveness  of  statistical  optimization  methods  in  achieving  restoration  and  super¬ 
resolution  depends  to  a  large  extent  on  the  probability  distribution  models  used  to  formulate  the 
optimization  criterion  (i.e.  the  probability  model  for  the  likelihood  function  or  for  the  a 
posteriori  distribution).  For  the  sake  of  obtaining  a  simple-to-implement  updating  algorithm,  one 
typically  employs  a  simple  distribution  function,  such  as  a  Poison  process,  to  model  these 
quantities.  Unfortunately  however,  processing  with  such  an  algorithm  can  yield  less  than 
optimal  performance  when  the  emission  process  underlying  the  particular  sensing  mechanism 
deviates  from  Poisson  statistics  and  requires  a  more  general  distribution  model  (such  as  a 
Rayleigh  distribution).  Use  of  more  general  distribution  functions  in  the  optimization  process 
will  result  in  corresponding  iterative  updating  schemes,  which  may  not  be  attractive  to 
implement  due  to  heavy  computational  demands.  Thus,  there  is  a  fundamental  conflict  between 
the  goals  of  implementation  simplicity  and  accuracy  of  representation  in  the  design  of  algorithms 
following  this  approach  forcing  one  to  sacrifice  one  or  the  other  of  the  two  objectives. 
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Yet  another  approach  for  image  restoration  that  is  gaining  greater  popularity  but  is 
fundamentally  very  distinct  from  the  statistical  optimization  methods  is  the  employment  of  a 
convex  set  theoretic  framework  in  which  several  constraint  sets  that  could  be  formulated  from 
known  a  priori  information  can  be  employed  during  iterative  processing  steps.  When  compared 
to  statistical  methods  for  image  restoration,  set-theoretic  methods  are  relatively  new.  The  earliest 
work  on  set-theoretic  estimation  in  a  general  context  can  be  traced  to  a  1 965  paper  by  Bergman 
[19],  The  Gerchberg-Papoulis  procedure  for  super-resolution  [4,5],  discussed  earlier,  can  be 
regarded  as  another  early  attempt  at  implementation  of  set-theoretic  estimation  to  signal 
reconstruction.  However,  it  is  the  introduction  of  a  new  iterative  procedure  for  an  organized 
synthesis  of  set  theoretic  estimates,  which  has  now  come  to  be  known  as  Projection  Onto  Convex 
Sets  (POCS),  that  marked  the  beginning  of  the  interest  in  this  approach  and  secured  its  steady 
growth.  Although  the  POCS  framework  was  first  developed  by  Gubin  et.  al.  [20],  and  further 
expanded  by  Lent  et.  al.  [21  ],  it  was  a  series  of  three  pioneering  papers  by  Youla  and  others  [22- 
24]  that  popularized  this  approach  among  the  scientific  and  engineering  community  and 
established  a  broad  conceptual  and  computational  basis  for  signal  recovery.  Some  useful 
properties  of  POCS  reconstructions  in  the  solution  of  deconvolution  problems  are  also  discussed 
by  Trussell  and  Civanlar  [25].  A  comprehensive  discussion  on  this  approach  together  with 
several  useful  extensions  can  be  found  in  a  tutorial  paper  by  Combettes  [26].  In  the  field  of 
image  processing,  more  interest  in  this  approach  has  so  far  been  shown  by  medical  imaging 
researchers  in  handling  such  problems  as  tomographic  image  reconstruction  [27,28]  and 
construction  of  3-D  animations  of  cardio-vascular  functions  [29],  etc. 

The  conflict  cited  above  that  exists  between  implementation  simplicity  (which  in  turn  is 
motivated  by  the  desire  to  execute  sufficient  number  of  algorithm  iterations  in  real  time  without 
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slowing  the  output  rate  of  sensor)  and  representation  accuracy  (which  is  motivated  by  the  desire 
to  utilize  as  much  of  the  relevant  and  accurate  information  as  feasible  during  the  processing 
steps)  has  encouraged  researchers  to  seek  new  methods  that  enable  one  to  handle  the  trade-offs 
more  satisfactorily.  In  this  report  we  shall  present  an  attractive  solution  towards  this  goal  by 
attempting  to  answer  the  question  -  “how  can  the  super-resolution  performance  of  a  simple-to- 
implement  algorithm,  such  as  the  Richardson-Lucy  iteration,  be  improved  by  use  of  additional 
knowledge  that  may  be  available?”  An  answer  to  this  question  is  suggested  by  the  inclusion  of 
intermediate  adjustment  steps  involving  implementation  of  convex  constraint  sets  that  model  the 
available  additional  knowledge.  Motivating  the  development  of  the  resulting  processing  schemes 
is  the  observation  that  set  theoretic  methods  offer  a  greater  flexibility  in  incorporating  available  a 
priori  information  into  the  restoration  process  than  what  one  can  achieve  by  attempting  to 
squeeze  all  of  this  information  into  a  single  probability  distribution  function  used  in  the  statistical 
optimization  process.  The  restoration  of  the  degraded  image  is  achieved  by  iteratively  projecting 
the  image  onto  the  various  convex  constraint  sets  that  model  the  different  types  of  a  priori 
information  one  may  have  in  the  restoration  problem  at  hand.  While  efforts  at  the  modeling  of 
constraint  sets  and  the  use  of  these  in  projection-based  set  theoretic  image  recovery  constitutes  a 
popular  direction  for  research  at  present,  it  seems  that  the  idea  of  combining  the  strong  points  of 
Bayesian  schemes  and  that  of  the  convex  set  based  methods  has  not  received  much  attention.  In 
fact,  development  of  processing  algorithms  that  intelligently  combine  the  strong  points  of 
Maximum  Likelihood  and  POCS  approaches  has  not  been  addressed  at  all.  A  primary 
contribution  of  this  project  is  the  design  of  hybrid  algorithms  that  attempt  to  leverage  the  strong 
points  of  both  ML  iteration  scheme  (simplicity  of  execution,  known  convergence  properties,  etc.) 
and  the  POCS  adjustments  (utilization  of  diverse  types  of  information  in  guiding  the  restoration 
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process,  faster  decay  of  restoration  error,  etc.).  The  development  of  specific  constraint  sets  that 
model  the  scene-derived  information  for  use  in  the  POCS  adjustment  steps  constitutes  yet 
another  important  contribution  of  this  work. 

1.5  Outline  of  the  Report 

The  organization  of  this  report  is  as  follows.  In  Section  2  we  discuss  some  basics  of  the  ML 
restoration  algorithm  and  describe  its  super-resolution  performance  when  applied  to  blurred 
images.  A  mathematical  analysis  of  the  dynamics  of  the  iterative  process  is  presented  to  give  an 
insight  into  the  properties  of  the  restored  image  as  well  as  to  gain  an  understanding  of  the 
possible  modifications  to  improve  the  super-resolution  performance.  In  Section  3  we  give  a  brief 
description  of  the  mathematical  background  for  set-theoretic  image  recovery.  We  also  propose 
two  novel  information  extraction  and  constraint  set  modeling  schemes,  and  their  application  in  a 
set  theoretic  image  recovery  framework.  In  Section  4  we  present  a  new  hybrid  approach  that  is 
tailored  to  combine  the  strong  points  of  a  Bayesian  scheme,  such  as  the  ML  algorithm,  and  those 
of  a  projection  based  set-theoretic  restoration  scheme.  Finally,  we  describe  the  performance  of 
two  new  super- resolution  algorithms  following  this  approach  by  presenting  results  of  a  few 
experiments  of  super-resolving  degraded  image  data.  In  order  to  emphasize  the  application  of 
these  algorithms  in  the  super-resolution  of  diffraction-limited  images  acquired  in  practice,  we 
will  specifically  include  results  from  experiments  of  processing  passive  millimeter-wave  images 
acquired  at  the  95  GHz  range1  from  state-of-the-art  millimeter- wave  radiometers.  Finally  in 
Section  6,  we  shall  present  a  summary  of  the  results  accomplished  during  this  project  and  offer 
some  concluding  remarks. 

1  It  is  to  be  noted  that  at  these  frequencies  the  signals  have  tremendous  penetrating  power  thus  providing  attractive 
imaging  capabilities  under  adverse  conditions  and  consequently  a  great  interest  exists  at  present  in  developing 
passive  sensors  operating  at  these  frequencies.  However,  due  to  Rayleigh  diffraction  limits,  the  resolution  in  the 
acquired  image  will  be  very  poor  offering  the  greatest  challenges  to  test  the  ability  of  restoration  and  super¬ 
resolution  algorithms. 
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2.  SUPER-RESOLUTION  BY  ITERATIVE  MAXIMIZATION  OF  LIKELIHOOD 
2.1  Restoration  and  Super-resolution  Performance  of  an  Iterative  ML  Algorithm 

As  noted  earlier,  iterative  image  restoration  algorithms  developed  using  a  maximum 
likelihood  (ML)  estimation  framework  have  attained  considerable  significance  in  recent  times 
for  their  super-resolution  capabilities.  The  basic  idea  underlying  these  methods  is  to  account  for 
the  statistical  behavior  of  the  emitted  radiation  at  the  level  of  individual  photon  events  by 
constructing  appropriate  object  radiance  distribution  models.  An  algorithm  that  affords  a 
convenient  digital  implementation  can  be  developed  starting  from  a  discretized  formulation  of 
the  model  for  the  imaging  process 

g(y)  =  ^O'.jO/M  +  hOO 

xeX 

where  /  (x)  denotes  the  object's  intensity  function  defined  on  a  region  X,  g(y)  denotes  the 
intensity  detected  in  the  image  defined  on  a  region  Y  and  h(x,  y )  denotes  the  point  spread 
function  (PSF)  of  the  imaging  sensor  and  n(y)  denotes  the  noise  present.  The  classical  image 
restoration  problem  is  to  find  the  object  intensity  estimate  f  (x)  given  the  data  g(y). 

There  exists  considerable  literature  on  dev  sloping  explicit  algorithms,  mainly  of  an  iterative 
nature,  for  handling  the  image  restoration  problem  within  a  statistical  framework  afforded  by 

such  a  formulation.  A  particularly  attractive  approach  is  to  obtain  a  ML  estimate  f(x),  i.e.  the 
object  intensity  estimate  that  most  likely  have  created  the  measured  data  g(y)  with  the  PSF 
process  h(y,x),  which  in  turn  is  developed  by  maximizing  an  appropriately  modeled  likelihood 
function.  Modeling  the  likelihood  function  is  basically  obtaining  a  goodness-of-fit  (GOF) 
quantity  for  the  measured  data,  since  the  likelihood  function  is  a  statistical  distribution  function 
p(g/J)  obtained  as  a  fit  to  the  relation  between  the  data  g(y)  and  the  object  f(x). 
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For  obtaining  a  digital  processing  algorithm,  assuming  that  the  image  to  be  processed 
consists  of  M  x  M  equally  spaced  gray  level  pixels  obtained  through  a  sampling  of  the  image 
field  at  a  rate  that  satisfies  the  Nyquist  criterion,  and  using  a  lexicographical  ordering  of  the 
signals  g,  f  and  n,  one  can  rewrite  (1)  as  a  convolution  of  two  one-dimensional  vectors 

h  =  [h(l),h(2),...,h(N)Y  and  /  =  [/(l),/(2),...,/(#)f  in  the  form 

g(i)  =  h(i )®  /(/)+  n(i),  i  =  1,2, ,N,  (2) 

where  N  =  M2  and  <8>  denotes  convolution.  Now  assuming  a  Poisson  distribution  for  p(g/f),  the 
estimate  of  the  object  f(x)  that  solves  the  optimization  problem 

/  =  arg  max  p[g/  f )  (3) 

can  be  developed  through  an  iterative  updating  scheme  given  by 

=  <4) 

where  fk  (j)  denotes  the  estimate  of  f(j)  constructed  at  the  k'h  iteration  and  <g>  denotes  discrete 

convolution.  Details  on  the  derivation  of  this  algorithm  can  be  found  in. the  literature  [4,5]. 

The  ML  algorithm  outlined  above  possesses  some  attractive  features  (such  as  guaranteed 
convergence,  non-negativity  of  estimates,  etc.)  that  are  useful  in  practical  implementation.  The 
optimization  framework  in  which  the  algorithm  is  developed  ensures  that  the  likelihood  function 
p(g/f)  monotonically  increases  over  successive  iterations  of  the  algorithm  and  hence  the 
processed  image  improves  in  quality  as  the  algorithm  proceeds. 

For  demonstrating  the  super-resolution  capabilities  of  this  algorithm,  some  representative 
results  of  applying  it  to  restore  blurred  images  is  shown  in  Figs.  1  and  2.  In  Fig.  la  is  shown  the 
blurred  version  of  a  simulated  object  comprising  of  a  series  of  concentric  disks  created  with  the 


19 


background  (dark)  at  intensity  value  zero  and  the  disks  (bright)  at  intensity  value  one.  For 
simulating  the  blurring  caused  by  a  diffraction  limited  imaging  sensor,  this  object  was  convolved 
with  the  PSF  of  a  low  pass  filter,  with  the  result  shown  in  Fig.  2a.  In  Fig.  2b  is  shown  the 
restored  version  of  this  image,  and  as  is  evident,  the  processing  removes  most  of  the  effects  of 
blurring  thus  attempting  to  restore  the  original  object.  To  demonstrate  the  extent  of 
super-resolution  achieved,  the  magnitude  spectra  of  the  blurred  image  and  the  restored  image  are 
shown  in  Figs.  3a  and  3b  respectively.  The  expanded  frequency  content  clearly  attests  to  the 
super-resolution  performance  of  the  algorithm. 


la  lb 


Figs.  2a  and  2b:  Blurred  “ Circles  ”  image  and  its  restoration  by  ML  algorithm 


2a  2b 


Figs.  3a  and  3b:  Magnitude  spectra  of  blurred  and  restored  Images 


To  illustrate  the  restoration  achievable  in  processing  practically  acquired  images  with  this 
algorithm,  the  results  of  an  experiment  are  presented.  Fig.  4a  shows  a  PMMW  image  (“Tank 
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Image”)  acquired  by  a  single  detector  radiometer  with  1  ft  diameter  aperture  operating  at  95 
GHz.  Figs.  4b  and  4c  respectively  show  the  restored  images  obtained  at  the  end  of  25  and  100 
iterations.  It  is  evident  that  the  resolution  has  considerably  improved  and  the  algorithm  is  quite 
powerful  in  producing  an  image  from  which  better  features  can  be  extracted  than  from  the 
original  for  any  further  information  exploitation  analysis,  such  as  target  classification.  For 
facilitating  real-time  implementation  capability  however,  one  desires  to  reduce  the  number  of 
iterations  needed  to  obtain  this  level  of  resolution  improvement,  or  to  achieve  even  better 
resolution  if  possible  with  less  computational  overhead.  As  will  be  shown  later  in  this  paper,  the 
POCS-assisted  ML  algorithm  developed  here  will  provide  this  capability. 


Figs.  4a,  4b  and  4c:  Acquired  PMMW  image  and  its  super-resolved  versions  after  25  and  100 

iterations  of  ML  algorithm 

2.2  Analysis  of  the  Dynamics  of  the  Iterative  Scheme 

Although  the  restoration  and  super-resolution  performance  of  the  algorithm  described  in 
Eq.  (4)  appear  quite  impressive  from  the  experimental  results  presented  above,  an  analysis  of  the 
dynamics  of  the  iterations  is  useful  to  obtain  a  greater  insight  into  the  restoration  performance.  In 
the  simulation  experiment  whose  results  are  shown  in  Figs.  2  and  3,  the  original  object  (“Circles” 
image  from  which  the  blurred  version  was  constructed)  is  known,  and  hence  a  plot  of  the 
“restoration  error”  (measured  as  some  appropriate  norm  of  the  deviation  between  the  original 
object  and  the  restored  image)  as  iterations  are  continued  can  provide  an  indication  of  the  speed 
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of  convergence  of  the  algorithm.  Also,  since  the  algorithm  is  developed  employing  a  maximum 
likelihood  framework,  at  each  successive  iteration  the  value  of  the  likelihood  function  can  be 


computed  from  noting  that  the  logarithm  of  the  likelihood  function  is  given  by 


f  N 


Uglf)=\ap(glf)  =  Y. 

k  =  I 

and  this  provides  an  alternate  way  of  displaying  the  dynamics  of  the  iterative  scheme.  Figs.  5a 
and  5b  respectively  show  the  behavior  of  the  Z,2  -norm  of  the  deviation  between  the  object  and 
the  processed  image,  and  the  increase  in  likelihood,  as  iterations  progress.  It  can  be  seen  that  the 
two  curves  display  an  inverse  behavior,  which  is  to  be  expected.  What  is  more  interesting, 
however,  is  to  observe  that  after  an  initial  steep  decay  of  the  restoration  error,  or  a  steep  increase 
in  the  likelihood,  the  changes  in  the  values  of  these  quantities  between  two  successive  iterations 
get  progressively  less  pronounced.  A  more  rigorous  mathematical  analysis  of  the  iterative 
scheme  is  useful  to  gain  further  insight  into  this  behavior. 


(5) 


Fig.  5  a.  Convergence  of  restoration  error  Fig.  5  b.  Likelihood  increase  with  progress  of 

iterations 
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The  ML  iterations  specified  by  Eq.  (4)  result  in  a  nonlinear  dynamical  process  of  the 


form 


LiU)  = 


(6) 


where  O (/*(/))  =  fkU) 


g(j ) 


.  The  convergence  of  the  iterations  can  hence 


be  examined  in  terms  of  the  equilibrium  points  of  this  dynamical  system,  or  equivalently  of  the 
“fixed  points”  fk  (j)  that  satisfy 

®(/» O'))  =  /»(/)•  (7) 

It  is  easy  to  see  that  this  condition  is  attained  at  the  £-th  iteration  when  the  restored  image 
satisfies  the  equality 

g(J)  =  fk(J)®h(J), 

since  the  quantity  inside  the  square  brackets  in  the  description  of  ®(/A(f))  becomes  1  (when  the 
PSF  is  symmetric  and  normalized,  which  in  turn  implies  that  1  <S>  h(J)  =  1 ).  Let  the  image 
estimate  vector  {/’(f)}  denote  this  converged  value.  Thus 

g(j)  =  fU)®Kj)  (8) 

specifies  the  fixed  points  j /*(y)}  of  the  ML  restoration  algorithm  given  by  Eq.  (2). 

A  careful  examination  of  the  condition  specified  by  Eq.  (8)  reveals  several  interesting 
facts.  First  of  all  observe  that  it  corresponds  to  the  condition  in  the  frequency  domain 

G(jco)  =  F\jco)H(jco)  (9) 

and  hence  is  satisfied  by  any  estimate  {/*(/))  that  restores  the  passband  spectrum.  Due  to  the 
band-limiting  operation  of  the  convolution  with  the  PSF  vector  ^e  image  estimate 
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{/*(/)}  that  satisfies  Eq.  (6)  is  not  unique.  Since  equality  outside  the  passband  of  the  sensor  is 
trivially  satisfied  (i.e.,  both  the  left  and  the  right  sides  of  Eq.  (9)  become  zero),  there  exist  an 
infinite  number  of  vectors  {/*(./)}  that  make  Eq.  (10)  hold.  These  vectors  are  however  related  in 
the  sense  that  each  of  them  has  a  specific  structure  that  permits  it  to  be  written  as  the  sum  of  a 
fixed  vector  {/^  (y)}  that  matches  the  object  spectrum  Fijon)  identically  below  the  diffraction 
limit  (i.e.,  for  co  <  coc)  and  another  vector  {/J>C/)}with  arbitrary  values  but  with  the  spectrum 

below  the  diffraction  limit  identically  zero.  In  other  words,  any  estimate  {/‘(y)}  that  can  be 
written  as 

r  u) =/;a)+/:  u)  m 

where  f'h  (  /)  and  f*h  (  /)  satisfy  the  frequency-domain  conditions: 

C i)H(jco)F'ph(jco)  =  G(jco ) 

(H)Kh  (ja>)  =  0  for  co  <  eoc,  (11) 

serves  as  a  fixed  point  of  the  ML  algorithm. 

It  is  now  simple  to  note  that  the  ML  algorithm  that  iteratively  strives  to  maximize  the 
likelihood  of  the  estimate  [  fk  (/)} ,  which  gives  rise  to  the  image  vector  |g(/)}  uPon  convolution 

with  the  PSF  vector  {/i(/)}  >  attempts  to  enforce  the  passband  restoration  condition  given  by  (/) 

above.  Hence,  when  an  estimate  {/*(/)}  that  satisfies  this  condition  is  obtained  at  the  &-th 

iteration,  the  algorithm  ceases  to  do  any  further  work.  Consequently,  convergence  of  the 
algorithm  is  to  be  interpreted  only  in  terms  of  restoration  of  the  passband  spectrum  of  the  image. 
In  other  words,  these  algorithms  principally  strive  to  restore  the  passband  as  accurately  as 
possible  with  the  following  resultant  effects: 
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(i)  Since  there  is  no  special  incentive  for  obtaining  the  part  {//*(/)}  that  matches  the  object 

spectrum  outside  the  diffraction  limit  specified  by  coc ,  this  vector  obtained  at  the  end  of  several 
iterations  (when  convergence  has  occurred)  could  have  elements  whose  contribution  to  the 
frequency  spectrum  is  much  smaller  than  the  contribution  of  the  other  part  .  In  other 

words,  the  “energy  in  the  frequency  extension”  could  be  very  small  (or  even  negligible). 

(ii)  When  passband  restoration  is  completed  (i.e.,  when  convergence  has  occurred),  no 

further  changes  in  the  estimate  {/*(/)}  are  likely,  since  the  likelihood  function  will  have 
attained  its  maximum  value. 

The  analysis  given  above  of  the  convergence  behavior  of  the  ML  iterations  is  very 
helpful  in  obtaining  an  insight  into  the  nature  of  the  multiple  fixed  points  of  the  algorithm  and 
the  effects  of  processing  digital  imagery  data.  Due  to  the  infinite  number  of  fixed  points  and  the 
corresponding  trajectories  that  describe  the  convergence  of  the  initial  estimates  to  these 
equilibrium  conditions,  obtaining  any  quantitative  estimates  of  the  convergence  rates  is  rather 
difficult.  The  problem  can  be  somewhat  simplified  by  examining  the  dynamics  of  the  “iteraiion 
error’,  or  change  in  the  estimate  during  successive  iteration  steps,  defined  by 

ek  O’)  =  fk  U)  -  /*-,  (/)  •  (12) 

The  behavior  of  this  error  is  governed  by  the  nonlinear  dynamical  process 

e.,,0)  -  L,U)-UJ)  -  <*»(/, (;))-/, O')  (13) 

which  can  be  re-written  as 

(w) 
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by  expressing  fk(j)  in  terms  of  ek(j) .  Eq.  (14)  represents  a  nonlinear  dynamical  system  whose 
trajectory  behavior  starting  from  an  initial  state  e0(j)  =  f0(J)  =  g(J)  is  of  interest.  The 
equilibrium  points  of  this  “error  system”  are  at  the  locations  where  T'(eA  (/))  =  0,  which  implies 
®(fkU))-fkU)  =  0>  and  hence  ek+](j)  =  0.  Thus  the  multiple  fixed  points  of  the  ML 
algorithm  {/'(y)}  (  on  the  N  -dimensional  space  for  the  estimate  f  )  are  mapped  into  the 

unique  equilibrium  point  { e(J )  =  o} ,  which  is  the  origin  of  the  //-dimensional  space  for  the 

iteration  error.  Hence  the  convergence  rate  for  the  restoration  can  be  studied  in  terms  of  the 
trajectory  behavior  of  the  “error  system”  given  by  Eq.  (14)  near  this  equilibrium  point. 

For  small  values  of  ek(j)  (i.e.,  near  convergence  to  this  equilibrium),  (/))  can  be 
linearized  by  a  Taylor  series  in  the  form 

xiJ(ek  (/))  =  Jek  (j)  +  [  Higher  Order  Terms  ]  (15) 

where  J  is  the  N  x  N  Jacobian  matrix  J  =  l  dek .  Ignoring  the  Higher  Order  Terms  in  the 
expansion,  the  iteration  error  dynamics  can  now  be  approximated  by  the  linear  system 

eMU)  =  Jek(J).  (16) 

Solving  this  system  of  equations  yields  the  trajectory  behavior  on  the  N  -dimensional  error 
space  which  can  be  described  by 

ekU)  =  Jke0U )  (17) 

where  Jk  is  the  state  transition  matrix  of  the  system  and  e0(j)  specifies  the  starting  error. 

From  the  convergent  nature  of  the  ML  algorithm  it  follows  that  all  eigenvalues  of  J  are 
inside  the  unit  circle,  i.e.,  |T,.(J)|  <  1  where  T;(J)  ,  i  =  1,2,  ...,  N ,  denote  the  eigenvalues  of  J. 
Consequently,  the  rate  of  decay  of  the  error  ek  ( j )  is  larger  during  the  first  few  iterations  of  the 
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algorithm  and  progressively  flattens  as  the  number  of  iterations  grow,  a  fact  that  was  confirmed 
by  the  simulation  experiments  earlier.  As  will  be  observed  later,  this  provides  the  motivation  to 
introduce  appropriate  adjustment  steps  after  executing  a  few  iterations  of  the  ML  algorithm  in 
order  to  re-gain  a  sharper  decay  of  the  error  by  re-playing  the  commencement  phase  of  the  ML 
algorithm  several  times. 
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3.  POCS  APPROACH  TO  RESTORATION  WITH  SCENE-DERIVED 

INFORMATION  CONSTRAINTS 

3.1  Basics  of  Convex  Set-theoretic  Restoration  of  Images 

An  approach  that  runs  parallel  to  but  is  very  distinct  from  the  statistical  optimization  method 
used  in  deriving  the  ML  algorithm  is  the  employment  of  a  convex  set  theoretic  framework  in 
which  several  constraint  sets  that  could  be  formulated  from  known  a  priori  information  could  be 
utilized.  In  order  to  highlight  the  differences  between  the  two  approaches,  starting  once  again 
with  the  discretized  model  of  the  imaging  process  given  by  Eq.  (1),  the  primary  objective  is  to 
model  each  piece  of  prior  information  VP/  that  is  known  as  a  closed  convex  constraint  set  Si  and 


seek  an  estimate  f  of  the  object  by  solving  the  optimization  problem 

f  =  arg  min7  £  W,J i  (/» s> ) 

/=! 


(18) 


where  L  denotes  the  number  of  constraints  defined,  wi  are  a  set  of  weights,  and  J ,  (/,  S, ) 

denotes  a  proximity  measure  that  determines  how  well  the  current  estimate  satisfies  the  i"‘ 
constraint.  The  proximity  measure  can  be  uniquely  defined  for  a  convex  constraint  scenario  as 
[26] 


MfA) 


minV.v, 


(19) 


where  f p  -  Pt  (/)  denotes  the  projection  of  /  onto  the  set  5,  .  It  is  evident  that  the  measure 

defined  above  quantifies  deviations  in  the  image  domain;  one  can  alternately  use  corresponding 
measures  that  quantify  deviations  in  the  frequency  domain  in  order  to  reflect  more  explicitly  the 
objectives  of  super-resolution. 

The  projection  operators  are  non-linear  mappings  with  several  analytical  properties  which 
can  be  exploited  to  develop  iterative  set-theoretic  algorithms  for  image  restoration.  Although 
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several  algorithms  were  developed  as  early  as  in  the  late  thirties,  the  potential  of  convex 
set-theoretic  methods  in  image  restoration  has  come  to  be  appreciated  since  the  advent  of  the 

method  of  projection  onto  convex  sets  (POCS)  [22-24], 

For  a  brief  description,  the  first  step  in  applying  the  method  of  POCS  to  an  image  recovery 
problem  is  to  define  a  closed  convex  set  for  each  of  the  a  priori  constraints  in  such  a  way  that  the 
members  of  the  set  are  consistent  with  the  associated  constraint  and  each  set  contains  the  actual 
image  distribution.  An  estimate  of  the  image  distribution  is  then  defined  to  be  any  member  of  the 
intersection  of  the  constraint  sets.  Finding  an  estimate  by  POCS  is  then  equivalent  to  the  problem 
of  finding  a  point  in  the  intersection  of  a  number  of  closed  convex  sets. 

Suppose  that  we  have  several  a  priori  known  constraints  that  can  be  associated  with  closed 
convex  sets  C, ,  i  =  1,  2,...,  L,  and  let  their  respective  projection  operators  be  denoted  Pi .  Then 

the  estimate  /„  generated  at  the  n"'  iteration  by  a  sequential  application  of  the  projections  on  the 

previous  iterate  ,  i.e. 

i,  =  r.P.->K- 2 . U,-<  (20) 

converges  to  an  estimate  /  in  the  intersection  set  C0=CinC2flC3 . DC.,.  It  is  further 

demonstrated  [26]  that  a  faster  convergence  to  an  estimate  /  lying  within  the  intersection  set 
can  be  achieved  if  instead  of  using  the  projection  operators  P,  in  the  construction  of  the  estimate, 
one  uses  “relaxed  projections”  7),  /  =1, 2, ....,  m  defined  by 

7]  =  (l  -  X)I  +  &P,  (21) 

where X  is  a  parameter  suitably  selected  within  the  range  0  <X<  2  and  the  estimate  /is 
constructed  iteratively  by  implementing  the  algorithm 


29 


f  -T  T  T 

Jn  1  tn1  m-\ x  /w-2 


(22) 


■Til-r 

For  evaluating  the  performance  of  these  algorithms  an  appropriate  metric  is  the  normalized 
proximity  function  given  by 

/;/ 

Zv,i/U) 

4>„  -  lOlogi! - .  (23) 

£»,w.a) 

1=1 

where  w(.  are  weights  selected  suitably  in  proportion  to  the  relative  importance  of  the  proximities 

to  individual  constraint  sets  the  final  estimate  is  desired  to  satisfy. 

The  projection  operations  executed  in  a  serial  manner  as  described  in  Eq.  (20)  (or  in  Eq. 
(22))  corresponds  to  enforcing  the  constraints  in  a  sequential  manner.  In  practice,  one  may  gain 
considerable  computational  time  benefits  by  applying  the  projections  in  parallel  and  constructing 
the  estimate  as  a  weighted  sum  of  the  individual  projections  as  will  be  described  later. 

3.2  Scene-derived  Information  and  Constraint  Set  Modeling 

Due  to  the  ill-posed  nature  of  the  super-resolution  problem,  it  is  essential  to  use  as  much 
prior  information  as  we  can  have  in  order  to  make  the  solution  more  and  more  well  conditioned. 
It  is  here  that  lies  the  advantage  of  projection  based  convex  set-theoretic  methods.  POCS  based 
image  recovery  schemes  are  extremely  flexible.  We  can  incorporate  any  a  priori  information  as 
long  as  there  exists  a  mathematically  tractable  projection  operator  for  evaluating  the  set 
estimates.  For  example  if  we  have  the  prior  knowledge  about  the  location  of  object  boundaries, 
then  we  can  model  it  as  a  convex  constraint  set  and  restore  the  object  boundary  edges.  Such 
restoration  of  the  object  boundaries  is  of  particular  interest  in  super-resolution  as  we  would  be 
restoring  the  high  frequencies  corresponding  to  the  object  boundary  edges. 
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Utilization  of  prior  information  for  convex  set-theoretic  super-resolution  of  images  can  be 
thought  of  as  a  two-step  process.  In  the  first  step  we  try  to  extract  appropriate  information  for  the 
image  to  be  restored.  The  second  step  is  primarily  involved  in  the  subsequent  modeling  of  that 
information  in  the  form  of  a  convex  set.  The  main  challenging  aspect  of  this  exercise  is  the 
extraction  of  reliable  information  from  given  images.  Constraint  sets  that  have  traditionally  been 
used  in  image  restoration  by  POCS  formulations  are  non-negativity  of  intensity  estimates,  known 
spatial  limits  on  objects  in  the  scene,  and  variance  of  noise  present  in  the  scene  being  estimated, 
all  of  which  are  fairly  easy  to  model  as  convex  sets  [28].  We  will  now  present  the  development 
of  two  new  constraint  sets  that  can  be  formulated  by  performing  simple  pre-processing 
operations  on  the  image  being  restored. 

Border  Constraint  Set: 

In  several  scenes  of  interest,  objects  of  interest  will  be  superimposed  as  foreground  against  a 
primarily  flat  background.  Thus,  a  sharp  transition  in  gray  level  intensity  exists  in  the  actual 
scene  imaged.  Due  to  blurring  caused  by  the  sensor,  this  sharp  edge  will  be  lost  and  image 
energy  from  foreground  region  flows  into  the  background  pixels,  thus  making  it  difficult  to 
decide  exactly  where  the  foreground  ends  and  the  background  begins.  Thus,  if  one  can  determine 
the  spatial  extent  of  the  object,  or  the  border  pixels  bounding  the  object,  one  can  sharpen  the 
edge  between  the  foreground  and  the  background.  In  the  spectral  domain,  this  sharpening 
corresponds  to  recovering  the  higher  frequency  components,  and  hence  super-resolution. 

Any  border  extraction  algorithm  based  on  edge  detection  processing  can  be  used  for 
developing  constraint  sets  that  permit  a  foreground-background  separation  of  the  image,  a 
popular  one  being  Canny’s  edge  detector  [30].  This  is  a  gradient-based  scheme  which  attempts  to 
design  the  edge  detector  by  employing  an  optimization  framework  in  order  to  achieve  the 


following  objectives:  (i)  maximize  the  signal-to-noise  ratio,  (ii)  achieve  good  localization  to 
accurately  mark  edges,  and  (iii)  minimize  the  number  of  responses  to  a  single  edge. 

A  version  of  this  scheme  that  is  particularly  suited  for  the  present  purposes  of  building  a 
constraint  set  for  POCS  restoration  will  now  be  outlined.  The  basic  idea  behind  the  procedure  is 
to  find  the  maximum  of  the  partial  derivative  of  the  image /in  a  direction  orthogonal  to  the  edge 
orientation.  To  execute  this  task,  one  can  create  a  two-dimensional  mask  for  the  edge  orientation 
by  convolving  a  linear  edge  detection  function  aligned  normal  to  the  edge  direction  with  a 
projection  function  aligned  parallel  to  the  edge  direction.  A  substantial  computational  savings 
results  if  both  functions  are  chosen  to  be  Gaussian  with  the  same  standard  deviation  cr  such  that 
the  mask  can  be  created  by  differentiating  normal  to  the  edge  orientation  as 


G 


n 


dG 
dn  ’ 


(24) 


where  G  is  the  two-dimensional  symmetric  Gaussian  function  G  -  exp 


.  Ideally,  n 


should  be  oriented  normal  to  the  edge  direction,  and  although  this  direction  is  not  known  a 
priori,  one  can  form  r  good  estimate  from  the  smoothed  gradient  direction 


A(G®  /) 
||A(G®/|' 


(25) 


where  <£>  denotes  convolution.  An  edge  point  is  then  defined  to  be  a  local  maximum  (in  the 

direction  of  ri)  of  the  operator  G„,  which  is  characterized  by  — G„  <8>  /  =  0 ,  or  equivalently 

dn 


A-(g®/)=  0. 

on 


(26) 
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Gaussian  smoothing  is  now  done  by  filtering  the  image  with  a  Gaussian  sliding  window,  with 
the  window  size  generally  chosen  to  be  6cr  with  the  two  dimensional  Gaussian  being  centrally 
located.  The  resulting  filtered  version  is  obtained  by  placing  this  mask  and  performing  a 
weighted  addition  of  the  pixels  which  lie  within  the  window,  with  the  weight  for  any  pixel 
chosen  as  the  value  of  the  Gaussian  at  that  location.  In  addition  to  smoothing,  this  step  also  helps 
in  reducing  the  number  of  spurious  edges  due  to  noise  that  may  be  present  in  the  image.  A 
judicious  selection  of  the  value  of  a  provides  the  desired  amount  of  smoothing,  with  larger 
values  providing  correspondingly  more  smoothing.  For  most  images,  we  found  that  a  value  of  cr 
selected  in  the  range  0.5  to  2  provides  a  good  choice. 

From  the  smoothed  image,  one  can  easily  compute  the  magnitude  and  the  direction  of  the 


gradient  as 


A/(U) 


df(i,  j)2  ,  8/(hjf 
dx  dy 


and 


©/(/,/)  =  arctan 


dx  dy 


(2?) 


using  the  approximations  ~  /(* +  ^  7)  ~  /0  ~  L  j)  an<^  ^~Qy~  f0*j  1)  • 

Once  the  gradient  values  are  obtained,  ordy  those  pixels  that  have  local  maxima  along  the 
direction  orthogonal  to  the  gradient  are  chosen  by  using  a  non-maximal  suppression  process 
which  involves  interpolating  the  gradient  field.  Edges  are  then  identified  as  those  pixels  that  have 
a  local  maximum  normal  to  the  gradient  direction.  Pixels  that  do  not  satisfy  this  criterion  are 
relegated  to  a  NO  EDGE  (0)  set,  while  those  that  pass  the  test  are  considered  as  POSSIBLE 
EDGE  (128)  candidates.  The  values  in  the  parentheses  correspond  to  gray  level  intensities. 

As  the  final  step  in  the  process,  the  possible  edges  are  linked  to  form  a  continuous  border 
by  implementing  a  double  threshold  linking  mechanism  which  involves  the  user  specifying  two 
threshold  parameters  thif,h  and  tUm .  Let  the  total  number  of  pixels  in  POSSIBLE  EDGE  be  Nedt.c . 
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Then  the  high  threshold  thiyh  corresponds  to  a  magnitude  level,  say  Mag high  such  that 
thi  h  *  Nali,c  number  of  pixels  have  a  magnitude  less  than  Maghjgh ,  which  can  be  determined 
from  a  magnitude  histogram  of  the  POSSIBLE  EDGE  pixels.  Once  the  high  threshold  Maghiyh  is 
determined,  the  low  threshold  is  then  set  to  be  Mag,im  =  thm  *  MagiUgh .  Any  pixel  in  POSSIBLE 
EDGE  that  has  a  gradient  magnitude  greater  than  Maghigh  is  classified  as  a  valid  EDGE  (0) 
point.  Any  pixels  connected  to  these  valid  EDGE  points  with  a  gradient  value  above  Maghw  are 

also  classified  as  EDGE  points.  In  other  words,  we  start  from  any  pixel  in  POSSIBLE  EDGE 
and  look  at  its  8-point  neighborhood  for  other  POSSIBLE  EDGE  pixels  that  satisfy  one  of  the 
two  criteria  stated  above.  At  the  conclusion  of  this  step,  the  algorithm  results  in  an  edge  map 
consisting  of  EDGE  (0)  pixels,  POSSIBLE  EDGE  (128)  pixels,  and  NO  EDGE  (0)  pixels. 
Selection  of  suitable  values  for  the  threshold  parameters  thigh  and  tlow  is  critical  for  a  continuous 
tracking  of  edges.  For  the  present  application  to  extract  border  pixels,  a  good  range  of  values  can 
be  determined  empirically  as  0.6  <  thi),h  <0.9  and  0.2  <  tbm.  <  0.4. 

The  border  tracking  performance  of  the  above  scheme .  can  be  illustrated  by  an 
experiment.  Fig.  6a  shows  the  input  image  used  for  analysis,  Fig.  6b  shows  the  edge  map 
obtained  and  Fig.  6c  shows  the  thresholded  edge  map  resulting  from  selection  of  threshold 
parameters  thigll  =  0.6  and  thm=  0.2. 
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6a 


6b 


6c 


Fig.  6.  Border  tracking  performance: 

6a.  Input  image;  6b.  Edge  map;  6c.  Thresholded  edge  map. 

Once  we  have  the  set  of  border  pixels  identified,  we  are  ready  to  proceed  with  the  second 
stage  of  the  border  information  modeling.  For  implementing  this  stage,  we  propose  a  method  of 
modeling  the  border  pixel  information  as  a  convex  constraint  set  and  subsequenty  define  a 
projection  operator. 

The  constraint  set  is  defined  as  one  that  enforces  a  background-foreground  separation  and 
can  be  mathematically  formulated  as 

Shurdcr  =  {/  €  E  :  foreground  of  image /  is  bounded  by  £  }  (28) 

where  £  denotes  the  set  of  border  pixels.  The  projection  operator  associated  with  this  set  exists 
provided  if  the  set  is  closed  and  convex,  which  can  be  easily  ensured  if  we  are  dealing  with  solid 
convex  shapes.  The  projection  operator  is  then  given  by 

■  ( i,j )  lies  inside  £ 

Pkank,  (/)  =  fp  where  fp  0»  j)  =  (  (29) 

0  :  (i,  j)  lies  outside  £. 

In  the  case  of  images  that  may  contain  multiple  objects,  a  challenge  in  the  application  of 
the  above  procedure  arises  in  making  a  decision  on  whether  a  pixel  lies  inside  one  of  the  borders. 
If  one  assumes  the  borders  to  be  convex  in  shape,  any  pixel  which  lies  inside  the  border  needs  to 
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be  completely  surrounded  by  it,  i.e.  if  we  move  up,  down,  left,  or  right  from  the  pixel  in 
question,  we  are  bound  to  hit  the  border  in  any  of  these  directions.  This  logic  may  lead  to  an 
erroneous  result  when  multiple  borders  are  present,  as  illustrated  in  the  scenario  depicted  in  Fig. 
6.  Hence,  not  only  do  we  need  to  have  intersection  with  the  set  of  border  pixels  as  we  move  in  all 
four  directions,  we  should  also  verify  whether  all  four  intersections  are  occurring  with  the  same 
set  of  border  pixels.  This  can  be  easily  handled  by  assigning  a  unique  label  (in  the  form  of  gray 
level  intensity)  to  each  of  the  different  sets  of  border  pixels.  As  illustrated  in  Fig.  7,  the  definition 
of  the  projection  operator  characterizes  pixel  at  location  P  as  a  valid  foreground  point  whereas 
both  Q  and  R  are  set  to  the  background  correctly. 


Fig.  7.  Projection  operation  for  enforcing  border  constraint  in  an  image  with  multiple  borders. 

To  obtain  a  feel  for  the  extent  of  restoration  and  spectrum  recovery  provided  by  this 
constraint  set,  a  controlled  experiment  was  conducted.  Fig.  8a  shows  a  simulated  image 
containing  multiple  objects  with  different  shapes  that  was  blurred  by  a  low  pass  filtering 
operation.  Fig.  8b  shows  the  blurred  image.  The  loss  of  spectral  content  due  to  blurring  can  be 
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seen  by  comparing  Figs.  8c  and  8d,  which  respectively  show  the  spectrum  of  the  original  image 
and  its  blurred  version.  The  border  map  obtained  from  an  application  of  the  procedure  described 
above  is  shown  in  Fig.  8e.  The  restored  image  obtained  from  a  single  enforcement  of  the  border 
constraint  is  shown  in  Fig.  8f  and  the  corresponding  spectrum  is  shown  in  Fig.  8g,  which  clearly 
illustrate  the  restoration  and  spectrum  recovery  achieved  from  this  constraint.  For  obtaining  a 


8e  8f  8g 


Fig.  8.  Restoration  and  super-resolution  of  a  blurred  image  using  border  constraint  set 
8a.  Test  image;  8b.  Blurred  image;  8c.  Spectrum  of  test  image;  8d.  Spectrum  of  blurred  image; 
8e.  Output  of  border  extraction  filter;  8f.  Restored  image;  8g.  Spectrum  of  restored  image. 

quantitative  evaluation  of  the  improvement  achieved,  the  L2  -norm  of  the  deviation  between  the 

original  image  and  the  restored  image  was  computed  as  18.69243,  which  compares  favorably 

with  the  corresponding  deviation  between  the  original  image  and  the  blurred  image  computed  as 

21.94160.  Achievement  of  this  level  of  restoration  from  a  single  application  of  the  constraint  is 

clearly  remarkable. 
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Template  Constraint  Set 

With  object  border  information  extracted  as  above,  we  cannot  restore  any  non-border  edges 
present  in  the  scene.  This  is  where  template-based  information  modeling  becomes  useful.  A 
requirement  for  utilizing  this  information  as  a  constraint  set  is  the  prior  knowledge  of  the 
presence  in  the  scene  of  a  certain  set  of  templates,  a  pattern  or  a  sub-image  that  may  occur  a 
multiple  number  of  times  in  the  scene  (a  right  angled  comer,  circular  wheels,  etc.).  A  template 
matching  filter  can  then  give  as  output  the  location  of  the  template  in  the  image. 

Several  procedures  for  design  of  template  matching  filters  are  discussed  in  literature.  Some 
early  methods,  which  are  quite  simple  to  implement,  are  based  on  image  subtraction  and 
correlation  matching.  In  one  of  these  methods,  the  absolute  difference  between  the  image  and  the 
template  is  used  to  determine  the  template  location.  While  computationally  less  demanding,  this 
method  assumes  perfect  knowledge  of  intensity  levels  and  works  only  in  a  stationary  imaging 
scenario.  Another  method,  which  is  more  robust  to  noise  and  image  variations,  attempts  to 
construct  the  filter  function  that  upon  convolution  with  the  image  containing  the  template 
generates  maximal  response  at  the  templat  locations.  The  main  drawback  in  this  method 
however  is  a  lack  of  proper  localization  of  the  template  occurrence  and  the  inability  to  account 
for  the  off-center  response  resulting  in  broad  peaks,  which  in  turn  can  interfere  with  accurate 
localization.  A  recent  and  more  efficient  procedure  for  designing  a  template  matching  filter 
involves  expansion  of  the  image  in  terms  of  nonorthogonal  basis  functions  [31],  which  we 
propose  to  employ  for  identifying  template  locations  and  construction  of  corresponding 
constraint  sets  for  restoration  processing.  Matching  by  expansion  is  based  on  expanding  the 
image  in  terms  of  basis  functions  which  are  all  self  similar  shifted  versions  of  the  template.  The 
procedure  is  briefly  outlined  in  the  following. 
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For  expanding  an  image  s(x,.y)of  size  NxN  in  terms  of  shifted  versions  of  the  template 
y) ,  one  uses  the  circulant  basis  functions 

^rc(x,y)  =  <f>(x-r,y-c),  r,c  =  0,1,....,jV-1  (30) 

in  terms  of  which  s(x,  y )  is  estimated  in  the  form 

Hx,  y ) = (*>  y)  (3i> 

In  (31),  Krc  denote  the  expansion  coefficients  which  are  to  be  found  to  minimize  the  deviation 
between  .?(x,.y)  and  .?(x,;y) .  By  a  lexicographic  ordering,  s(x,y) ,  i(x,  j/) ,  and  y/K(x,y)  can  be 
expressed  as  N 2  -dimensional  vectors  V,V,  and  T,  respectively  in  order  to  formulate  the 
criterion  for  finding  the  coefficients  Krc  as  one  of  minimizing  the  squared  deviation  given  by  the 
inner  product  <(V -V),(V -V)>.  Now  utilizing  the  orthogonality  principle  that  the 
approximation  error  V  -  V  is  orthogonal  to  the  basis  functions,  one  obtains  the  N 2  equations 
<(V- h%>  =  0,  /  =  0,1,2,...., jV2  - 1 .  (32) 

For  facilitating  the  goal  of  evaluating  the  expansion  coefficients  Kfc,  these  equations  can  be 

rewritten  more  compactly  in  a  matrix  form  as 

=  W  (33) 

where  denotes  the  N2  x  2 autocorrelation  matrix  whose  (/,  j)  -th  element  is  ^'J//,VF/^  ,VF  = 
j  J  and  Kis  a  matrix  consisting  of  the  expansion  coefficients Krc .  If  the  basis 
functions  selected  are  linearly  independent,  is  positive  definite  and  Krc  can  be  uniquely 

determined  by  solving  Eq.  (33).  Solution  of  Eq.  (33)  by  inverting  is  computationally 
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prohibitive  in  most  cases,  and  hence  Ben-Arie  and  Rao  [31,32]  suggested  an  alternative  simpler 
method  obtained  by  treating  the  expansion  problem  as  an  equivalent  deconvolution  problem. 

For  the  discrete  image  s(x,y)  of  size  NxN  to  be  expanded  by  a  set  of  self-similar  basis 


functions  y/ij(x,y)as 

s(x>y) = (34) 

if  the  basis  functions  i//{j{x,y)  can  be  densely  arranged  in  discrete  equidistant  locations,  we  have 


Vij  (x,y)=  </>(x,  y)®8{x-i,y-  j)  =  </>{x-i,y-  j ) 


(35) 


and  hence  (34)  can  be  re-written  as  the  discrete  convolution 


s(x,  y)  =  _ Ky<f>(x  -i,y-  j )  =  k{x,  y)  0  </>(x, y).  (36) 

The  problem  of  finding  the  expansion  coefficients  kh  can  now  be  handled  as  reversing  the 

convolution  operation  in  Eq.  (36).  A  particularly  simple  method  is  to  design  a  Wiener  filter  for 
the  case  when  the  image  s(*,.y)  is  assumed  to  contain  a  noise  signal  A(x,y) ,  i.e. 

s(x,y)  ==  K{x,y)  0  </>(x,y)  +  Z(x,y) .  (37) 

The  Wiener  filter  is  then  given  by 


0(©,,®2) 


SJ^oy 2) 


S„  (fl), ,  co 2 )  =  0(<y,  ,0) 2 )  *  Scc  (©, ,  co2 )  ;  S„  (©, ,  a>2 )  =  5CC  (©, ,  <u2  )||0(<y, ,  <y2  )||2  +  Su  (<y, ,  ©2 ).  (3  8) 


In  Eq.  (38),  S^  and  5^  respectively  denote  the  spectral  densities  of  the  signal  5(x,_y)and 
noise  A(x,  y) ,  and  denotes  the  cross  spectral  density  of  5(x,y)  and  the  coefficients  rc(x,y). 
O (<y,,rti2)  is  the  Fourier  transform  of  the  template  <p(x,y)  and  the  symbol  *  denotes 
conjugation.  In  an  ideal  scenario,  the  discrete  function  /r(jc?  _y)  will  ha\e  a  value  1  for 
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x  =  x0 ,  y  =  y0  and  will  be  zero  elsewhere.  Hence  the  spectral  density  of  x(x,y )  should  be  a 
constant,  say  Kv  Moreover,  if  the  additive  noise  is  assumed  to  be  white  gaussian,  its  power 
spectral  density  is  also  a  constant,  say  K2.  Hence,  Eq.  (38)  can  be  rewritten  as 


©(<»,,  ©2  )  = 


\$>{co„co2f+K 


(39) 


where  K  =  ^~ 

K, 


9a 


Fig.  9.  Template  location  in  a  test  image 

9a.  Template;  9b.  Test  image  containing  template;  9c.  Plot  of  expansion  coefficients 


For  an  illustration  of  the  application  of  this  filter  to  identify  template  locations,  Fig.  9a 
shows  a  template  that  occurs  at  two  locations  (10,10)  and  (150,130)  in  the  test  image  shown  in 
Fig.  9b.  Fig.  9c  shows  a  plot  of  the  magnitudes  of  the  expansion  coefficients  in  various  columns, 
which  correctly  depict  sharp  peaks  at  the  two  relevant  columns  where  the  template  is  present. 

Once  the  locations  of  a  template  in  the  image  are  obtained,  one  can  proceed  to  design  a 
constraint  set  and  a  corresponding  projection  operator  for  use  in  the  POCS  restoration  scheme. 
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For  a  template  of  size  mxm  identified  at  the  ( k,l )  pixel  location  by  the  filter,  the 

constraint  set  can  be  defined  as 

Slempl!lll,  =  {/  e  E :  f(k  +  i,l  + /)  =  ^(/, ;) ,  where  /, y  e  [0,1,2 . m  - 1]  }  (40) 

(k,l)  e  {  Chc :  Set  of  coordinate  locations  of  the  template  occurrences  } . 

The  projection  operator  for  this  set  is  given  by 

<t>{Uj)  :  Uj  e  [0,l,2,...,wi-l]  } 

(/)  =  fP  where  fP(k  +  i,l  +  j)  =  {  (41) 

/(/,/ )  :  elsewhere. 

The  template  matching  filter  as  described  above  performs  a  binary  decision  regarding  the 
existence  or  non-existence  of  a  specific  template  at  a  given  location  in  the  image,  and  outputs  a 
peak  value  at  the  template  location.  If  multiple  templates  are  present  in  a  given  image,  a  bank  of 
filters,  each  constructed  for  a  specific  template,  can  be  employed  in  parallel  for  obtaining  the 
template  locations.  Alternately,  it  is  possible  to  modify  the  filter  design  such  that  peaks  of 
different  magnitudes  corresponding  to  different  templates  are  obtained  in  the  filter  output  thus 
eliminating  the  need  for  multiple  filters.  Details  of  these  modifications  are  omitted  here  and  can 
be  found  in  [31,32j. 

3.3  Restoration  of  Blurred  Images  using  Scene-derived  Constraint  sets 

In  this  section  we  shall  briefly  describe  the  restoration  and  super-resolution  performance  of 
two  set-theoretic  estimation  algorithms  using  the  constraint  sets  developed  in  the  last  section  by 
employing  a  POCS  framework.  In  addition  to  the  two  sets  Shonk,r  and  S,  lale  described  by  Eqs. 

(28)  and  (40),  we  will  use  the  following  three  constraint  sets: 

(i)  SHiaK  =  {/  €  S  :  255  >  /(/J)  >  0  ,i,j  =  1,2,.  ...N  }  (42) 

(ii)  Scxlem  =  {/  e  E  :  f(i,j )  =  0  if  (i,j)  £  £  }  (43) 
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and 


(iii)  S„„.  “  (/  e  2 :  =  G(W).(*./) «  f  )•  <44> 

While  the  constraint  imposed  by  SIIU,K  is  self-explanatory,  SuxM„  imposes  a  spatial  extent 
constraint  based  on  any  prior  knowledge  of  a  support  £  for  the  actual  image  distribution,  and 
5  imposes  a  constraint  in  the  frequency  domain  by  setting  the  Fourier  transform  F(k,l)  of 
the  image /equal  to  the  Fourier  transform  G(k,l)  of  the  original  image  g  within  the  passband  of 
the  sensor  £ .  It  may  be  noted  that  these  constraint  sets  are  the  ones  that  have  been  traditionally 
employed  in  POCS-based  image  reconstruction  work  in  the  past  [22-24],  The  corresponding 
projection  operators  can  be  defined  as: 

if  0< /(/,/><  255 

(0  (f)  =  fp  where /,(/,»=  {  0  ,  if  0  (45) 

255  ,  if  255. 

/(/,;),  if  (/,/)  e  £ 

(ii)  P,r,a„  (/)  =  fP  where  fP  (h  j)  =  I  (46) 

0  ,  if 

G(k,l),  if  (k,l)  €  £ 

(iii)  Pieman  (/)  =  fp  where  FP  (k>‘ 0  =  (  ^ 

F(k,l) ,  if  (k,l)e£. 

These  three  operators  together  with  the  operators  Phnnk.r  and  Plaiiplale  defined  in  Eqs.  (29)  and  (41) 
form  the  set  of  operators  that  will  be  used  in  a  POCS  framework  for  restoration  and  super¬ 
resolution. 

For  designing  the  specific  algorithm  that  employs  these  constraints  one  has  several  options. 
One  option  is  to  use  a  sequential  implementation  of  the  projections  Pi  at  each  iteration  step  as 
described  by  Eq.  (20)  in  order  to  construct  the  next  iterate.  Alternately,  as  noted  earlier,  one  may 
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use  the  relaxed  projections  Tj ,  constructed  from  Pi  as  described  in  Eq.  (21),  in  a  sequential 

manner  to  implement  the  algorithm  described  in  Eq.  (22).  Improved  results  together  with  faster 
implementation  of  the  algorithm  can  be  obtained  from  a  parallel  execution  of  the  projections  at 
each  iteration.  We  shall  hence  discuss  this  method  in  a  greater  detail  in  the  following. 

As  noted  before,  in  the  traditional  implementation  of  the  projections  one  enforces  the 
constraints  in  a  sequential  manner,  which  in  general  can  result  in  large  computation  times, 
particularly  if  a  number  of  constraint  sets  are  used.  A  more  efficient  implementation  will  involve 
employing  a  parallel  processing  architecture  in  which  the  projections  onto  the  different  constraint 
sets  at  each  iteration  can  be  simultaneously  executed  and  the  next  iterate  is  computed  as  a 
weighted  sum  of  these  projections.  One  version  of  this  approach,  popularly  referred  to  as  the 
Method  of  Parallel  Projections  (MOPP)  [33,34],  has  been  found  in  our  investigations  to  possess 
attractive  convergence  properties.  In  this  scheme,  an  elementary  iteration  consists  in  projecting 
the  current  estimate  simultaneously  onto  selected  sets  and  forming  a  relaxed  convex  combination 
of  the  projections. 

For  a  brief  description  of  the  method,  given  the  image  to  be  restored,  denoted  as  the  starting 
image  /(0)  for  commencing  the  POCS  iterations,  and  two  numbers  e  s  (0,l)  and  M  a  positive 
integer,  MOPP  implements  the  recursion 


together  with  the  following  conditions: 

(i)  X„  e  [e,2-e]  and  (ii)  =  1,  where  0  *  /„  c  I  and  I  c  UA/j0lf„+* , 

iel„ 

I  denoting  the  set  of  constraints  used. 
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For  demonstrating  the  restoration  performance  of  the  MOPP  algorithm,  we  shall  present  the 
results  of  a  simulation  experiment  in  Figs.  lOa-b.  The  blurred  image  to  which  the  method  was 
applied  is  shown  in  Fig.  10a  and  its  restored  version  after  20  iterations  of  MOPP  algorithm  using 
the  border  and  template  constraint  sets  is  shown  in  Fig.  10b.  The  normalized  proximity  measure 
given  by  Eq.  (8)  was  used  as  the  performance  measure. 


10a  10b 


Fig  1 0.  Restoration  of  a  blurred  image  with  MOPP  algorithm 
10a.  Blurred  image;  10b.  Restored  image 
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4  A  HYBRID  SUPER-RESOLUTION  ALGORITHM  COMBINING 
POCS  AND  ML  APPROACHES 


It  is  evident  from  the  discussion  in  the  last  section  that  the  POCS  approach  provides  an 
attractive  alternative  to  the  statistical  optimization  approach  for  restoring  degraded  images.  In 
particular,  for  super-resolving  poorly  acquired  images  by  the  POCS  approach  the  flexibility 
offered  for  constructing  constraint  sets  and  utilizing  them  during  processing  is  particularly 
noteworthy.  It  is  to  be  noted  that  the  types  of  information  used  in  the  restoration,  viz.  border 
constraints  and  template  constraints,  are  not  possible  to  utilize  during  implementation  of 
algorithms  derived  following  statistical  optimization-based  approaches,  for  instance  during 
execution  of  ML  super-resolution  algorithm  discussed  in  Section  2.  However,  statistical 
optimization  methods  in  general,  and  ML  algorithm  in  particular,  have  a  number  of  advantages 
(as  discussed  before)  which  one  may  not  like  to  sacrifice.  The  strong  points  of  these  two 
complementary  approaches  can  be  combined  to  produce  hybrid  schemes  that  offer  very  attractive 
restoration  and  super-resolution  performance.  Two  such  hybrid  algorithms  will  be  presented  in 
this  section. 

4.1  POCS-assisted  ML  algorithm 

In  this  scheme,  we  implement  the  ML  iterations  given  by  Eq.  (2)  as  the  main  image 
restoration  scheme  and  include  a  POCS  adjustment  of  the  ML  estimate  for  enforcing 
constraints  after  every  cycle  of  ML  iterations  is  completed.  More  specifically,  we  execute  L 
iterations  of  POCS  adjustment  after  each  cycle  of  K  iterations  of  the  ML  algorithm.  Thus, 
each  iteration  cycle  of  the  combined  algorithm  applies  an  ML  iteration  cycle  comprising  K 
number  of  ML  iterations  followed  by  a  POCS  adjustment  cycle  comprising  L  number  of 
POCS  iterations.  The  combined  algorithm  is  run  for  a  total  of  N  iterations.  It  is  easy  to  see  that 
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this  algorithm  in  effect  performs  a  constrained  maximization  of  the  likelihood  function  in  an 
iterative  manner.  The  implementation  details  are  summarized  in  the  following  table. 


J.  Read  the  blurred  image  data  (g). 

A 

2.  Set  initial  conditions  to  commence  ML  iterations,  f0  =  g. 

3.  Perform  Ml  restoration  implementing  the  updating  algorithm  given  in  Eq.  (2). 

4.  Compute  'he  likelihood  and  l2  -norm  values  at  the  end  of  each  iteration. 

5.  Repeat  steps  3-4  for  K  iterations. 

6.  Perform  POCS  adjustment  using  MOPP  scheme  and  implementing  updating  algorithm  given  in  Eq.  (48). 

7.  Repeat  step  6  for  L  iterations. 

8.  Repeat  steps  3-6  for  N  iterations. 

9.  Save  processed  image  and  compute  likelihood  and  l2  -norm  values. 


Table  1.  Simulation  Flow  of  POCS-assisted  ML  Algorithm 
For  evaluating  the  performance  of  this  new  algorithm,  we  conducted  a  number  of  simulation 
experiments.  The  results  of  one  experiment  will  be  briefly  described.  In  this  experiment,  the 
POCS-assisted  ML  algorithm  was  used  to  restore  the  blurred  image  shown  in  Fig.  10a  by 
employing  the  same  two  constraints  as  discussed  in  the  last  section  (  viz.,  Border  Constraint  and 
Template  Constraint)  for  POCS  adjustment.  The  algorithm  was  implemented  with  values  of  K=  5 
and  N=3  ( K  being  the  number  of  ML  iterations  executed  before  a  POCS  adjustment  and  N  being 
the  total  number  of  iterations  for  the  combined  cycle).  However,  we  varied  the  value  of  L,  the 
number  of  POCS  adjustment  iterations  within  each  cycle  from  two  to  ten  and  examined  the 
effect  it  has  on  the  performance  of  the  ML  algorithm.  Figs.  1  la  and  lib  respectively  show  plots 
of  the  variation  of  the  likelihood  function  and  the  l2  -norm  of  restoration  error  computed  in  the 
spatial  domain  against  the  total  number  of  ML  iterations  executed  with  the  value  of  L  varied 
from  L  =  0  (pure  ML  algorithm  with  no  POCS  adjustments)  to  L  =  2,  6,  and  10.  It  is  evident  that 
the  hybrid  scheme  outperforms  the  ML  algorithm  without  any  POCS  adjustments.  An  item  of 
special  interest  to  note  is  the  jump  in  the  likelihood  function  as  the  algorithm  switches  from  ML 
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iteration  cycle  to  the  POCS  adjustment  cycle.  A  corresponding  fall  in  the  l2-norm  of  the 
restoration  error  is  also  observed  at  these  points.  Furthermore,  the  amount  of  increase  in  the 
likelihood  (and  the  decrease  in  /,  -norm  of  the  restoration  error)  improves  with  the  number  of 
iterations  of  POCS  ( i.e.  value  of  L )  executed.  However,  it  is  also  interesting  to  see  that  most  of 
the  improvement  due  to  the  inclusion  of  POCS  adjustments  comes  during  the  first  few  iterations 
of  POCS  and  beyond  that  there  is  not  much  change  in  the  performance. 


Fig.  lib.  Plot  of  l2  -norm  of  Restoration  Error 
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4.2  ML-assisted  POCS  algorithm 

An  alternate  way  of  creating  a  hybrid  restoration  algorithm  is  to  use  POCS  updating  as  the 
main  restoration  scheme  and  employ  execution  of  ML  updating  for  the  intermediate 
adjustment  steps.  We  shall  refer  to  this  algorithm  as  ML-assisted  POCS  algorithm  in  further 
discussion.  More  specifically,  we  execute  K  iterations  of  ML  adjustment  after  each  cycle  of  L 
iterations  of  the  POCS  adjustment  implemented  using  the  MOPP  approach  as  in  Eq.  (48). 
Thus,  each  iteration  cycle  of  the  combined  algorithm  now  applies  a  POCS  iteration  cycle 
comprising  L  number  of  POCS  iterations  (implementing  the  updating  in  Eq.  (48))  followed  by 
a  ML  adjustment  cycle  comprising  K  number  of  ML  iterations  (implementing  the  updating  in 
Eq.  (2)).  The  combined  algorithm  is  run  for  a  total  of  N  iterations.  The  implementation  details 
are  summarized  in  the  following  table. 

/.  Read  the  blurred  image  data  (g). 

2.  Set  initial  conditions  to  commence  POCS  iterations,  f0  =  g. 

3.  Perform  POCS  restoration  using  MOPP  framework  and  adaptive  relaxation  implementing  the  updating 
algorithm  given  in  Eq.  (48). 

4.  Compute  the  likelihood  and  l2  - norm  values  at  the  end  of  each  iteration. 

5.  Repeat  steps  3-4  for  L  iterations. 

6.  Perform  ML  adjustment  implementing  updating  algorithm  given  in  Eq.  (2). 

7.  Repeat  step  6  for  K  iterations. 

8.  Repeat  steps  3-6  for  N  iterations. 

9.  Save  processed  image  and  compute  likelihood  and  l2  - norm  values. 

Table  2 .  Simulation  Flow  of  ML-assisted  POCS  Algorithm 

The  performance  of  this  algorithm  was  also  evaluated  by  applying  it  to  the  same  restoration 
problem  discussed  earlier  of  restoring  the  blurred  image  shown  in  Fig.  10a.  The  POCS  updating 
was  implemented  by  employing  the  same  two  constraints  as  discussed  in  the  last  section  (  viz., 
Border  Constraint  and  Template  Constraint).  The  algorithm  was  implemented  with  values  of  L=  5 
and  N=3  (L  being  the  number  of  POCS  iterations  executed  before  a  ML  adjustment  and  N  being 
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the  total  number  of  iterations  for  the  combined  cycle).  In  this  case,  the  value  of  K,  the  number  of 
ML  adjustment  iterations  within  each  cycle,  was  varied  to  examine  the  effect  it  has  on  the 
performance  of  the  POCS  restoration.  Figs.  12a  and  12b  respectively  show  plots  of  the  variation 
of  the  normalized  proximity  measure  (computed  as  in  Eq.  (23))  and  of  the  l2-norm  of  restoration 
error  computed  in  the  spatial  domain  against  the  total  number  of  POCS  iterations  executed  with 
the  value  of  K  varied  from  K  =  0  (pure  POCS  algorithm  with  no  ML  adjustments)  to  K  =  2,  and 
5.  Once  again  it  is  evident  that  the  hybrid  scheme  outperforms  the  POCS  iteration  scheme 
without  the  ML  adjustments.  Of  particular  interest  is  to  note  the  downward  break  in  the  two  plots 
as  the  algorithm  switches  from  the  POCS  iteration  cycle  to  the  ML  adjustment  cycle. 
Furthermore,  the  amount  of  performance  improvement  increases  with  the  number  of  iterations  of 
ML  adjustment  (i.  e.  value  of  K)  executed  between  two  POCS  cycles. 


Fig.  12a.  Plot  of  Proximity  Measure 
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Fig.  12b.  Plot  of  l2-norm  of  Restoration  Error 


Hoe  does  this  algorithm  compare  with  the  other  hybrid  algorithm  presented  earlier?  To 
answer  this  question,  an  examination  of  the  /2-norm  plots  in  Figs,  lib  and  12b  can  be  made 
since  the  same  blurred  image  was  processed  by  the  two  algorithms  using  the  same  sets  of 
parameters.  At  the  end  of  15  iterations  of  the  main  restoration  algorithm,  the  ML-assisted  POCS 
yielded  a  smaller  value  of  the  restoration  error  than  the  POCS-assisted  ML  algorithm.  The  latter 
algorithm  started  with  a  larger  initial  error  but  the  POCS  adjustment  step  caused  a  sharper  decay 
in  the  error  compared  to  that  produced  by  the  ML  adjustment  step  in  the  ML-assisted  POCS 
algorithm.  One  should  also  take  into  consideration  that  the  ML-assisted  POCS  algorithm 
required  more  computation  time  for  the  same  number  of  iterations  than  the  other  scheme.  Thus  in 
practical  implementations  where  processing  time  (and  hence  the  number  of  frames  that  can  be 
processed  within  a  given  time)  is  of  particular  concern,  one  may  prefer  the  POCS-assisted  ML 
algorithm  over  the  other  hybrid  scheme.  On  the  other  hand,  when  quality  of  the  final  restored 
image  is  of  importance,  the  ML-assisted  POCS  algorithm  seems  to  be  a  better  choice. 
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4.3  Application  of  POCS-Assisted  ML  Algorithm  to  Restoration  of  PMMW  Images 


An  application  of  POCS-assisted  ML  algorithm  has  been  made  to  the  super-resolution  of 
PMMW  images  with  significant  performance  improvements  noted  in  comparison  to  the 
application  of  ML  algorithm.  For  providing  an  illustration  of  this,  we  show  in  Fig.  13a  a  PMMW 
image,  popularly  referred  to  as  “Pattern  Image”,  collected  by  a  95  GHz  camera  with  a  3  in 
lens/horn  built  at  the  Army  Research  Laboratory  (ARL)  at  Adelphi,  MD.  This  image,  of  size 
65x33  pixels,  has  over  the  years  become  a  sort  of  benchmark  image  for  testing  the  performance 
of  different  super-resolution  algorithms.  More  details  on  this  data  can  be  found  in  various  ARL 
publications,  one  of  which  is  the  recent  paper  by  Silverstein  [18]. 

Fig.  13b  shows  the  restored  image  after  15  iterations  of  the  ML  algorithm.  The  application  of 
the  POCS-assisted  ML  algorithm  to  the  same  image  with  5  ML  iterations  and  3  POCS 
adjustments  alternately  implemented  for  2  cycles  (i.e.  K  =  5,  L  =  3,  and  N  =  2)  resulted  in  a 
restored  image  shown  in  Fig.  13c.  Comparing  the  images  in  Figs.  13b  and  13c,  one  may  note  the 
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Fig.  13.  Results  of  processing  PMMW  “pattern”  image  by  ML  and  POCS-assisted  ML 

algorithms. 

13a.  Acquired  Image;  13b.  ML  Processed  Image;  13c.  Image  processed  by  POCS-assisted  ML 
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superior  edge  restoration  and  the  absence  of  processing  artifacts,  which  clearly  illustrate  the 
strengths  of  the  present  hybrid  algorithm.  It  is  to  be  noted  that  these  performance  gains  have 
resulted  in  addition  to  a  reduction  in  the  total  number  of  ML  iterations  (only  10  in  POCS-assisted 
ML  compared  to  15  in  ML),  which  is  the  more  computationally  expensive  part  of  the  super¬ 
resolution  process. 

As  a  further  comparison  of  the  super-resolution  performance,  the  spatial  frequency 
distribution  plots  were  also  compared.  Fig.  14  shows  the  variation  of  the  normalized  intensity  (in 
dB)  plotted  as  a  function  of  the  frequency  for  a  cut  through  the  spectrum  of  the  acquired  image 
(shown  in  dotted  lines)  and  also  the  corresponding  variations  for  the  ML  processed  image 
(shown  in  dark  solid  line)  and  for  the  POCS-assisted  ML  processed  image  (shown  in  fainter  solid 
line).  It  is  easy  to  see  that  while  the  ML  algorithm  has  recovered  some  frequencies  beyond  the 
passband  of  the  sensor,  the  recovery  of  spectral  components  by  the  POCS-assisted  ML  algorithm 
is  significantly  more  extensive. 

As  a  final  illustration  of  the  performance  of  the  POCS-assisted  ML  algorithm  in  super¬ 
resolving  PMMW  images,  we  show  in  Fig.  15  below  the  restoration  of  the  “Tank”  image 
discussed  earlier  in  Section  2.  This  is  a  resultant  of  processing  with  the  POCS-assisted  ML 
algorithm  using  5  ML  iterations  alternating  with  5  POCS  adjustment  iterations  implemented  for 
3  cycles,  which  corresponds  to  a  computation  time  roughly  equivalent  to  that  required  to  execute 
15  ML  iterations.  Comparing  the  restoration  obtained  here  with  that  shown  in  Fig.  4c,  which  is 
the  restored  image  after  100  ML  iterations,  one  can  clearly  appreciate  the  benefits  resulting  from 
the  hybrid  approach  presented  in  this  work. 
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Magnitude  in  db 


Spectral  Data(ARL  "Pattern"  Benchmark) 


Fig.  14.  Comparison  of  spatial frequency  spectra  of processed  and  unprocessed  images 
(Dotted  line  -  original  blurred  image;  Dark  solid  line  -  restoration  by  ML  algorithm; 
Fainter  solid  line  -  restoration  by  POCS-assisted  ML  algorithm) 


Fig.  15.  Restoration  of  “Tank”  image  by  POCS-assisted  ML  algorithm 


5.  CONCLUSIONS 


The  major  contributions  from  this  project  are  the  development  of  scene-based  information 
sets  for  use  as  constraints  within  a  set-theoretic  estimation  framework  for  image  restoration  and 
the  design  of  two  new  hybrid  algorithms,  called  POCS-assisted  ML  algorithm  and  ML-assisted 
POCS  algorithm.  These  algorithms,  which  iteratively  apply  a  cycle  of  ML  estimation  iterations 
followed  by  a  cycle  of  POCS  adjustment  iterations,  combine  the  strong  points  of  the  two 
approaches  and  hence  possess  a  number  of  implementation  benefits.  For  facilitating  the  POCS 
adjustment,  two  novel  information  sets  that  may  be  derived  from  the  image  being  processed  were 
also  developed.  These  sets,  together  with  the  traditionally  applied  constraints,  such  as  non¬ 
negativity  of  intensity,  spatial  extent  of  objects,  and  noise  variance,  provide  a  rich  set  of 
information  constraints  that  can  be  advantageously  utilized  to  speed  up  the  convergence  of  the 
ML  estimation  process.  Several  experimental  results  were  presented  to  describe  the  performance 
of  the  new  algorithms  in  super-resolving  PMMW  images.  The  superior  restoration  of  the  object 
features  observed  in  the  image  domain  and  the  significant  extrapolation  of  spatial  frequencies 
observed  in  the  spectral  domain  lead  us  to  conclude  that  the  POCS-assisted  ML  algorithm  and 
the  ML-assisted  POCS  algorithm  are  perhaps  the  most  powerful  super-resolution  algorithms 
available  today  for  restoration  of  imagery  acquired  from  diffraction-limited  sensing  operations. 
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